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We study quarkonium correlators and spectral functions at zero and finite temperature using the 
anisotropic Fermilab lattice formulation with anisotropy £ = 2 and 4. To control cutoff effects we 
use several different lattice spacings. The spectral functions were extracted from lattice correlators 
with Maximum Entropy Method based on a new algorithm. We find evidence for the survival of 
fS quarkonium states in the deconfined medium till relatively high temperatures as well as for 
dissolution of fP quarkonium states right above the deconfinement temperature. 
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INTRODUCTION 



The study of heavy quarkonium correlators at zero and 
finite temperature is interesting for several reasons. First, 
quarkonia are hadrons whose properties and structures 
are best understood; to the first approximation their 
properties can be described in terms of non-relativistic 
potential. Also due to their small sizes they can provide a 
bridge between perturbative and non-perturbative QCD. 
On general grounds it is expected that quarkonia will 
melt at temperatures somewhat higher than the decon- 
finement temperature as a result of modification of inter- 
quark forces (color screening). Thus their properties at 
finite temperature can serve as an indicator of in-medium 
modification of inter-quark forces and help to under- 
stand the phenomenon of non-Abelian Debye screening in 
Quark Gluon Plasma. Furthermore, it was suggested by 
Matsui and Satz [l| that color screening at high tempera- 
ture will lead to quarkonium suppression which can serve 
as a signal of Quark Gluon Plasma formation in heavy 
ion collisions. Quarkonium properties at finite temper- 
ature have been extensively studied in potential models 
[3, i i 1 i, i 0, S ©, M, HO, El, M, ■ However, 
the validity of potential models at finite temperature is 
doubtful, ft is more appropriate to study modifications 
of quarkonium properties at finite temperature in terms 
of spectral functions. First principle calculations of the 
quarkonium spectral functions at finite t emp erature have 
become available only recently [HI, [H, Il7| . In the cal- 
culations of charmonium spectral functions which have 
been performed so far, the systematic errors were not 
well controlled. Although some features of the spectral 
functions extracted so far have been assigned to be due 
to lattice discretization effects, it has not been examined 
in detail which features of the charmonium spectral func- 



tions are physical and which are merely artifacts of the 
finite lattice spacing and analysis methods. 

The aim of the present paper is twofold. First, we 
would like to study the properties of the quarkonium 
spectral functions in a large range of lattice spacings to 
control systematic effects both due to finite lattice spac- 
ing and limited number of data points. Second, we want 
to study quarkonium correlators at non-zero tempera- 
ture to see what happens to different charmonium states 
above deconfinement. Having at hand results from differ- 
ent lattice spacings allows to check the cutoff dependence 
of these results. In addition to charmonium correlators 
we also calculate bottomonium correlators and spectral 
functions. Some preliminary results of this study have 
been presented in [H, El HE H3 • 

The rest of the paper is organized as follows. In sec- 
tion II we discuss general properties of quarkonium cor- 
relators at finite temperature. In section III we describe 
the analysis tools for the reconstruction of the quarko- 
nium spectral functions. In section IV we discuss our 
lattice setup used in the calculation of the charmonium 
correlators. In section V we discuss charmonium spectral 
functions at zero temperature, emphasizing the system- 
atic uncertainties of the analysis. Section VI deals with 
the temperature dependence of charmonium correlators. 
In section VII we discuss the problem of reconstruction 
of charmonium spectral functions at finite temperature. 
Section VIII contains some results on charmonium corre- 
lators at non-zero spatial momenta. In section IX and X 
we discuss our results for bottomonium correlators and 
spectral functions. Finally section XI contains our con- 
clusions. The reader not interested in technical details 
may skip sections II, III and IV. 
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II. MESON CORRELATORS AND SPECTRAL 
FUNCTIONS 

In this section we discuss the relation between the Eu- 
clidean meson correlators and spectral functions at finite 
temperature. It is straightforward to take the zero tem- 
perature limit. 

Most dynamic properties of the finite temperature sys- 
tem are incorporated in the spectral function. The spec- 
tral function (7h{potP) f° r a given mesonic channel H in 
a system at temperature T can be defined through the 
Fourier transform of the real time two point functions 
Z? > and D < or equivalently as the imaginary part of the 
Fourier transformed retarded correlation function 1231 , 
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In the present paper we study local meson operators 
of the form 



with 



J H (t, x) = q{t, x)T H q(t, x) 



T# = 1,75,7^ 757^,7^ 



(3) 



(4) 



for scalar, pseudo-scalar, vector, axial-vector and tensor 
channels. The relation of these quantum number chan- 
nels to different meson states is given in Tab. U 

The correlators D^ < \xq, x) satisfy the well-known 
Kubo-Martin-Schwinger (KMS) condition [22j 



D^(x ,x) = D < (x + i/T,x). 



(5) 



Inserting a complete set of states and using Eq. ([5]) , one 
gets the expansion 



(n\J H {0)\m)\H\p^-kl + k™) 



(6) 



where Z is the partition function, and k n ^ m ^ refers to the 
four-momenta of the state \n(m)). 

A stable mesonic state contributes a S function-like 
peak to the spectral function: 
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where ttih is the mass of the state. For a quasi-particle 
in the medium one gets a smeared peak, with the width 
being the thermal width. As one increases the temper- 
ature the width increases and at sufficiently high tem- 
peratures, the contribution from the meson state in the 



spectral function may be sufficiently broad so that it is 
not very meaningful to speak of it as a well defined state 
any more. The spectral function as defined in Eq. (|6|) 
can be directly accessible by high energy heavy ion exper- 
iments. For example, the spectral function for the vector 
current is directly related to the differential thermal cross 
section for the production of dilepton pairs [23| : 
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Then presence or absence of a bound state in the spectral 
function will manifest itself in the peak structure of the 
differential dilepton rate. 

In finite temperature lattice calculations, one calcu- 
lates Euclidean time propagators, usually projected to a 
given spatial momentum: 

G H (r,p) = J d 3 xe l P s (T T J H (T,x)J H (0,6)) (9) 

This quantity is an analytical continuation of D > (xo,p) 

G H (T,p)=D>(-iT,p). (10) 

Using this equation and the KMS condition one can easily 
show that Gh(t, p) is related to the spectral function, Eq. 
0]), by an integral equation (see e.g. appendix B of Ref. 

my- 



G(r,p) 
K(lo,t) 



dioa(u;,p)K(uj, t) 



cosh(w(r- X/2T)) 



sinh(w/2T) 



(11) 



This equation is the basic equation for extracting the 
spectral function from meson correlators. Methods to do 
this will be discussed in the next section. Equation (JTTJ) 
is valid in the continuum. Formally the same spectral 
representation can be written for the Euclidean correlator 
calculated on the lattice G lat (T,p). The corresponding 
spectral function, however, will be distorted by the effect 
of the finite lattice spacing. These distortions have been 
calculated in the free theory [U HE] • 



III. BAYESIAN ANALYSIS OF MESON 
CORRELATORS 

The obvious difficulty in the reconstruction of the spec- 
tral function from Eq. (|TT|) is the fact that the Euclidean 
correlator is calculated only at O(10) data points on the 
lattice, while for a reasonable discretization of the in- 
tegral in Eq. (fTT|) we need 0(100) degrees of freedom. 
The problem can be solved using Bayesian analysis of the 
correlator, where one looks for a spectral function which 
maximizes the conditional probability P[a\DH] of hav- 
ing the spectral function a given the data D and some 
prior knowledge H (for reviews see [26|, H3]). Different 
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TABLE I: Meson states in different channeis 



Bayesian methods differ in the choice of the prior knowl- 
edge. One version of this analysis which is extensively 
used in the literature is the Maximum Entropy Method 
(MEM) [Hill. It has been used to study different cor- 
relation functions in Quantum Field Theory at zero and 
finite temperature [StSEl HE HE M, M, HI, H HE 
HE HE HE HE EE SlElP In this method the basic 
prior knowledge is the positivity of the spectral function 
and the prior knowledge is given by the Shannon - Janes 
entropy 

cr(ui) — m(oj) — a(u>) ln( — 7— r) 

The real function m{oj) is called the default model and 
parametrizes all additional prior knowledge about the 
spectral functions, e.g . such as the asymptotic behav- 
ior at high energy [26l [30j . For this case the conditional 
probability becomes 

1 



S = du> 



P[a\DH] = exp(--x 2 +aS), 



(12) 



with x 2 being the standard likelihood function and a 
a real parameter. In the existing MEM analysis of the 
meson spectral functions the Bryan's algorithm was used 

Here we will introduce a new algorithm described be- 
low. To maximize (TT2|) we have to solve the following 
equation 



Sa(uj) 
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-X -aS 
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(13) 



If we perform N measurements at points t%, . . . ,tn with 
results Gi (i = 1 . . . N) , and the data model Gi [a] is 
described by (fTl"|) . then we obtain 



N 



In 



m(u) 



= 0, (14) 



where Cy is the correlation matrix. With the notation 



1 N 

3=1 



Gi 



we can write this equation as 
o(uj) — m (w) exp 



' N 



(15) 



Since Si itself depends on a(u>), this is just a re- 
parametrization of the original problem. However, by 
substituting back this form into (fT4|) we obtain 



N 

E 

i=l 



N 



3=1 



0. (16) 



Since the functions K(u),Ti) are linearly independent, 
this equation can hold only if the expressions in the 
square brackets are zero. Multiplying by the correlation 
matrix it reads 



N 

3=1 



(17) 



Since the relation a[s] is known one can solve these equa- 
tions. 

We can go even a step further by recognizing that 



G j [o-] = 
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This allows us to define a function 
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a N 

— SiCijSj 

i,3'=l 
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Using this function p?|) can be written as 
dU 



ds. 



= 0. 



(18) 



(19) 



(20) 



It can also be shown that d 2 U / dsidsj is positive defi- 
nite. Therefore we reformulated the MEM problem to 
the task of minimizing U . This is a minimization prob- 
lem in number-of-data dimensions, which is not a more 
difficult task than applying a x 2 method. So we should 
expect that our MEM-code runs as fast as a \ 2 code, as 
it is indeed the case. We used the Levenberg-Marquant 
method to perform the minimization which in most cases 
is stable enough to find the minimum. 
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It is worth to make connection between our method 
and the Bryan algorithm [2{|. In both cases the true 
problem is number-of-data dimensional - in more dimen- 
sions the problem would be under-determined. To find 
the relevant subspace, the Bryan algorithm uses singular 
value decomposition, while wc find the same relevant sub- 
space by exact mathematical transformations. Although 
the method of identifying the subspace is different, the 
result is the same, and in both cases one proceeds with 
solving the original problem in this restricted subspace. 
The advantage of the new algorithm is that it is more sta- 
ble numerically when we reconstruct quarkonium spec- 
tral functions at zero temperature. We will discuss this 
issue in more detail in section Ivl 

To demonstrate the reconstruction power of MEM and 
also to reveal the role of noise in the data we give some 
illustrative examples. We use the simple uncorrelated 
noise model 



5G{t) = 6rG(r), 



(21) 



and vary b. For the given spectral function we gen- 
erate mock data for N T — 32 time slices. The terms 
proportional to a introduce the dependence on the prior 
knowledge, so it is advantageous to choose for it small- 
est value possible; the typical value in the analysis was 

1 -8_ 10 -12 

The first is a typical continuum spectral function con- 
sisting of a Dirac delta and a continuum part. In Fig. [T] 
the original spectral function can be seen together with 
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FIG. 1: The continuum-like spectral function and its MEM- 
reconstruction in presence and in absence of noise, respec- 
tively. 

the MEM-reconstruction in the presence and absence 
of noise, respectively (the magnitude of the noise was 
b = 1CP 4 ). The matching is not perfect, but all the main 
characteristics are present. The ability of MEM to re- 
construct finite width spectral functions has been demon- 
strated in a similar analysis with the Breit-Wigner ansatz 
used to model the ground state instead of the Dirac delta 
[2lj . The simple noise model in this analysis assumes a 



diagonal correlation matrix. In real situations the corre- 
lation matrix is not diagonal. We will discuss this case 
later in the course of the analysis of the lattice correla- 
tors. 



IV. LATTICE FORMULATIONS AND 
PARAMETERS FOR CHARMONIUM PHYSICS 

To study the charmonium system at zero and finite 
temperature we use the anisotropic Fcrmilab formulation 
described in Ref. [43[ . In our study we use the quenched 
approximation and use the standard Wilson action in the 
gauge sector for which the relation between the bare £o 
and the renormalized anisotropy £ = a s /at is known in a 
wide range of the gauge coupling j3 = 6/g 2 [44|. For the 
heavy quarks we use the anisotropic clover action [43| 



rn + v t p 



Wilson , YjL 



2 ( ^sw £ atsFts + ~F L £ ass,Fss ' ) 



Here the Dirac operator is defined as 
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with 



Un{x)ip{x + n) - Ufa - n)ip(x - n) 



A^{x) = 

U^(x)ip(x +fi) + Ufa- fj,)ip(x - pi) - 2ip(x) 



(24) 



(25) 



Furthermore, a^^, = {7^,7^} and the field strength ten- 
sor is defined as 

f»Ax) = -I[Q»»-QU (26) 

±QA X ) = U^x)U,(x + fi)Ufa + v)Ufa) + 

u v (x)ufa - a + Wfa - ®u M (x - A) + 

Ufa - (i)Ul{x - fi - !>){7 Al (x - fc— v)U v (x — v) + 
Ul{x - i^U^x - v)U v {x +fx- i>)Ufa). 

(27) 

As in Ref. [43| we fix v s = 1 and use tadpole improved 
values for the clover coefficients. The remaining param- 
eters of the action, the bare heavy quark mass Too and 
the bare velocity of light v t are fixed non-perturbatively. 
The bare quark mass is fixed by requiring that either 
the pseudo-scalar (?y c ) mass or the spin averaged IS mass 
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are equal to the experimental value. To fix v t we study 
the dispersion relation E(j>)o{ r\ c and require that the 
velocity of light defined by [43j 



^/l6E 2 { Pl ) - f5£ 2 (;po) - E 2 ( P2 ) 



Pn 



f 12Ap 
2nn/N a ,Ap = 2tt/N s 



(28) 



is equal to one. Here N s is the spatial size of the lattice. 
For crosscheck we also calculate the velocity of light given 
by the simpler definition 



y/£ 2 (pi)-£; 2 (Po) 
Ap 



(29) 



Both definitions give consistent results. In our study we 
use £ = 2, 4 and (3 = 5.7, 5.9, 6.1, 6.5 corresponding to 
temporal lattice spacings a^ 1 — 1.905 — 14.12 GeV. The 
bare anisotropy £o f° r a given renormalized anisotropy 
was calculated in Ref. [4J| for (3 — 5.5 — 6.5. To set 
the scale for the lattice spacing we use the tradition phe- 
nomenological value ro = 0.5 fm for the Sommer scale 
[45l ]. The Sommer scale ro has also been calculated for 
anisotropic Wilson action for f3 = 5.5 — 6.1 [46]. For 
larger j3 values we use extrapolation based on RG in- 
spired (3 function [47] 
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R(6) 
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with R(f3) = (6o/3)" w/(2b « ) exp(-l/(26 )/3) being the 
standard 2-loop beta function. The extrapolation is done 
for (3 = 6.5 using (3 = 5.9, 6.0 and 6.1. To crosscheck the 
results we also did the extrapolation including the point 
at f3 = 5.8 and the difference gives the estimate of the 
error on the extrapolated ro. Alternatively one can es- 
timate the lattice spacing from the difference between 
the mass of 1 P\ state and the spin averaged IS 1 mass: 
AMI^Pi — 15). To a very good approximation this mass 
difference is not affected by fine and hyperfine splitting 
and thus is not very sensitive to quenching errors. It 
was found that close to the continuum limit the lattice 
spacing determined from AM( l P\ — IS) is different from 
that determined from r by 10% [H, |48| if we use the 
phenomenological value r — 0.5 fm. Using the value of 
r = 0.469(7) determined in full QCD [4§] would give a 
value for AAf( 1 Pi — IS*) splitting which is closer to the 
experimental one, however, the AM(2S — IS) splitting 
would be even further away from the experimental value 
[48j |. This problem is due to the quenched approxima- 
tion. The parameters of simulations are summarized in 
Tab. HQ For £ = 2 they are identical to those of Ref. [Z^]. 
For £ = 4 the values of the clover coefficients have been 
taken from Ref. (5lj |. and other two parameters, mo and 
v t were fixed non-perturbatively as described above. 

By studying the large distance fall-off of the correlators 
we determine the mass of the ground state in each quan- 
tum number channel. The pseudo-scalar ground state 
masses from the fit are also given in Tab. [TTJ 



V. CHARMONIUM SPECTRAL FUNCTIONS 
AT ZERO TEMPERATURE 

In this section we discuss the calculations of the char- 
monium spectral functions at zero temperature for dif- 
ferent channels. The continuum meson current in Eq. [3l 
Jh is related to the lattice current as 



J H = Z H a s ipT H ip, 



(31) 



where ip is the lattice quark field in Eq. ([22)1 . The 
renormalization constant Zh can be calculated in per- 
turbation theory or non-perturbatively. Unfortunately 
for anisotropic clover action they have not been calcu- 
lated. Since in this paper we are mainly interested in 
controlling systematic errors in the spectral functions and 
the temperature dependence of the charmonium correla- 
tors the exact values of the renormalization constants 
are not important. Nonetheless, when we compare the 
spectral functions calculated at different lattice spacing 
in the plots it is convenient to take into account the ef- 
fect of Zh- We do this in an ad- hoc way by choosing 
Zh such that the continuum part of the spectral func- 
tion is roughly equal to its free value at large energies. 
The values of Zh are given in the Appendix. In what 
follows we will discuss charmonium spectral functions at 
zero spatial momentum p = 0. 



A. Pseudo-scalar and vector spectral functions at 
zero temperature 

In this subsection we discuss our results on charmo- 
nium spectral functions obtained using MEM. To recon- 
struct the spectral functions we used the algorithm de- 
scribed in section [TTTl We also compare the results with 
the Bryan algorithm for the pseudo-scalar spectral func- 
tion for (3 = 6.1, £ = 4. The problem with the Bryan 
algorithm is that it does not work well for charmonium 
correlators if the time extent is sufficiently large, which 
is the case at low temperatures; the iterative procedure 
does not always converge. For instance at j3 = 6.1, £ = 4 
and 16 3 x 96 lattice we could get the spectral functions us- 
ing the Bryan algorithm only when using t max = 24 data 
points in the time direction. With the new algorithm 
there is no restriction on t max which can be as large as 
N-t/2. The comparison of the spectral functions obtained 
with the two algorithm is given in the Appendix. 

When we analyze the correlation function using MEM 
we need two inputs : the default model m(w) and the 
parameter a. Thus the spectral function we get depends 
both on a and the default model. We investigated the 
a dependence of the spectral functions. For sufficiently 
large a the spectral function is very smooth and shows 
no peaks. As we decrease a the spectral function shows 
more and more peaks. These features are also shown in 
the Appendix. It has been shown that one can construct 
the probability P[a|Z?m] of ha ying some a for given data 
and default model m(u) [2^, [2i|. Typically, for high 



6 



p 


5.7 


5.9 


6.1 


6.1 


6.5 


N 2 S x N t 


8 3 x 64 


16 3 x 64 


16 3 x 64 


16 3 x 96 


24 2 x 32 x 160 


(£> &) 


(2,1.6547) (2,1.6907) (2,1.7183) (4,3.2108) 


(4,3.31655) 


ro/a a 


2.414(8) 


3.690(11) 


5.207(29) 


5.189(21) 


8.96(4) 




1.905 


2.913 


4.110 


8.181 


1/119 


^ sw 


2.138 


1.889 


1.7614 


1.9463 


1.7054 


^ sw 


1.3252 


1.2055 


1.1431 


1.0984 




mo 


0.51 


0.195 


0.05 


0.05 


-0.025 




1.01 


1.09 


1.12 


1.25 


1 1 Q 

i.iy 


M(?7c)[GeV] 


3.029(1) 


3.052(1) 


2.994(2) 


2.983(3) 


3.031(9) 


c(0) 


1.000(2 ) 


0.984(3) 


0.984(3) 


1.002(4) 


1.001(7) 


L a [fm] 


1.66 


2.17 


1.54 


1.54 


1.34 


configs 


2000 


1560 


930 


500 


630 



TABLE II: Simulation parameters for charmonium at zero temperature. Also shown here is the mass of the r\ c state and the 
renormalized speed of light c. 

statistics data this probability is a smooth function of 
a and has a maximum at some a = a max (26l . l3lT | . To 
eliminate the a dependence the spectral functions calcu- 
lated at fixed a were averaged over alpha with the weight 
P[a\Dm] 0, [MESS HI. In the present analysis 
we simply calculate the spectral function at a max - In 
the new algorithm used in this paper the sub-space, in 
which the search for the spectral function (i.e. the maxi- 
mization of the conditional probability P[a\DH] ) is per- 
formed, is chosen by selecting different data on the cor- 
relator G(t),t = 0, 1,2, ...,N t /2 which go into the anal- 
ysis. Although the dimension of this sub-space can be 
as large as N t /2 in practice it is limited by statistics. In 
our analysis the dimension of this subspace was tuned to 
the largest possible value giving a smooth a-dependence 
of P[a|-Dm] with a unique maximum. We use different 
type of default models, all of them being smooth function 
of uj, m(u>) = mQio 2 , m(u>) — h = const as well as the 
form given by the free spectral function calculated on the 
lattice [24j . To reduce the sensitivity to the lattice ar- 
tifacts at short time separation arising from anisotropy, 
in the reconstructed spectral function we do not include 
the data on the correlator for t < £ in our analysis. This 
was also done in Ref. It turns out, however, that 

including these points does not alter significantly the re- 
sults. 

The zero temperature spectral functions for three dif- 
ferent lattice spacings are summarized in Fig. [5] Here we 
used the simple default model m(oj) = 1 (in this paper 
we always give the default model in units of the spatial 
lattice spacing ) . To get a feeling for the statistical errors 
in the spectral functions we calculate its mean value in 
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FIG. 2: The pseudo-scalar spectral function at zero temper- 
ature for three finest lattice spacings. The horizontal line 
corresponds to the spectral function in the free massless limit 
at zero lattice spacing. 



some interval /: 



a = 



Jj du>a(u)) 



(32) 



Then we calculated the error on a using standard jack- 
knife method. These errors are shown in Fig. [21 where 
the length of the intervals are shown as horizontal error 
bars. As one can see from the figure, the r) c (lS) can be 
identified very well. The second peak is likely to cor- 
respond to excited states. Because of the heavy quark 
mass the splitting between different radial excitations is 
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FIG. 3: The default model dependence of the pseudo-scalar 
spectral function at the finest lattice spacing (/3 = 6.5). In 
the insertions the default models corresponding to different 
spectral functions are shown. 



small and MEM cannot resolve different excitations in- 
dividually but rather produces a second broad peak to 
which all radial excitation contribute. This can be seen 
from the fact that the amplitude of the second peak (i.e. 
the area under the peak) is more than two times larger 
than the first one. Physical considerations tell us that 
it should be smaller than the first amplitude if it was a 
2S state. When comparing amplitudes and peak posi- 
tions from MEM analysis and from double exponential 
fits we find very good agreement for the first peak and a 
fair agreement for the second peak. The details of this 
comparison are given in the appendix. This gives us con- 
fidence that at zero temperature charmonium properties 
can be reproduced well with MEM. 

For energies larger than 5 GeV we probably see a con- 
tinuum in the spectral functions which is distorted by 
finite lattice spacing. In particular the spectral function 
is zero above some energy which scales roughly as aj 1 . 
Note that for ui < 5GeV the spectral function does not 
depend on the lattice spacing. 

One should control how the result depends on the de- 
fault model. In Fig. [3] we show the spectral function for 
three different default models. One can see that the de- 
fault model dependence is significant only for u) > 5 GeV. 
This is not surprising as there are very few time slices 
which are sensitive to the spectral functions at lo > 5 
GeV, while the first peak is well determined by the large 
distance behavior of the correlator. 

We would like to know how well we can reconstruct the 
spectral function using MEM. To answer this question 
we consider a model spectral function consisting of four 
peaks and the perturbative continuum. We used semi- 
realistic masses and amplitudes calculated from potential 
model [ill] . From this spectral function we generated 
mock data for the correlator G mock (i) at lattice spacing 
aT 1 = 14.12GeV and N t = 160. The correlation matrix 
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FIG. 4: The comparison of spectral functions obtained from 
MEM with the input spectral function (thick black lines). 



used in the analysis of this mock data was defined as 



C 



mock 



G(i)G(j) 



G 



n>< ck /j\Qmock 



(33) 



Here Cij and G(i) are the correlation matrix and me- 
son correlators calculated in the pseudo-scalar channel 
on 24 2 x 32 x 160 lattice at = 6.5. 

The result of the MEM analysis of this mock data are 
shown in Fig. [4J where the four delta functions and the 
continuum of the model spectral function are shown as 
thick black lines. The first peak is reproduced quite well. 
The three radial excitations show up as one broad peak. 
Indeed, the area under the second peak is 1.05GeV 3 while 
the input was 1.20GeV 3 . The continuum is only repro- 
duced correctly when the default model has exactly the 
same form as the continuum at large u>. 

So far we only discussed the pseudo-scalar channel. We 
also calculated the spectral function in the vector channel 
defined as 



cry (uj) 



(34) 



The zero temperature vector spectral functions are shown 
in Fig. [5] for the three finest lattice spacings. The con- 
clusions which can be derived from this figure are similar 
to the ones discussed above for the pseudo-scalar chan- 
nel. The first peak corresponds to the J/i/)(lS) state, 
the second peak most likely is a combination of 2S and 
higher excited states, finally there is a continuum above 5 
GeV which is, however, distorted by lattice artifacts. The 
similarity between the pseudo-scalar and vector channel 
is, of course, expected. The lower lying states in theses 
channels differ only by small hyperfine splitting. 



B. Spectral functions for P states 

We calculated the spectral functions in the scalar, 
axial- vector and tensor channels which have the IP char- 
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FIG. 5: The vector spectral function at zero temperature for 
three finest lattice spacings. The horizontal line is the spectral 
function in the free massless limit. 
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FIG. 6: The scalar spectral function at zero temperature for 
three finest lattice spacings. 
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FIG. 7: The axial-vector (top) and tensor (bottom) spectral 
functions at zero temperature for three lattice spacings. 



spectral function is again strongly distorted. 



VI. CHARMONIUM CORRELATORS AT 
FINITE TEMPERATURE 



monia as the ground state. The scalar spectral functions 
reconstructed using MEM are shown in Fig. [6] The 
first peak corresponds to \c0 state, but it is not resolved 
as well as the ground state in the pseudo-scalar channel. 
This is because the scalar correlator is considerably more 
noisy than the pseudo-scalar or vector correlator. For the 
two finest lattice spacings there is a second peak which 
may correspond to a combination of excited P states. 
Above u> > 5 GeV we see a continuum which is strongly 
distorted by lattice artifacts and probably also by MEM. 
In the Appendix we show the comparison of the ampli- 
tudes and the masses of the ground state from MEM and 
two exponential fits. 

The spectral functions in the axial-vector and tensor 
channels are shown in Fig. \7\ They look similar to the 
scalar spectral functions. As in the scalar channel the 
first peak is less pronounced than in the case of S wave 
charmonium spectral functions, and it corresponds to Xci 
and h c state, respectively. The continuum part of the 



As we have seen in the previous section the reconstruc- 
tion of the spectral functions from lattice correlators is 
difficult already at zero temperatures. At finite temper- 
ature it is even more difficult to control the systematic 
errors in the spectral functions reconstructed from MEM. 
This is because with increasing temperature the maximal 
time extent r max is decreasing as 1/T. Also the number 
of data points available for the analysis becomes smaller. 
Therefore we are looking for a method which can give 
some information about the change of the spectral func- 
tions as the temperature is increasing. The temperature 
dependence of the spectral function will manifest itself 
in the temperature dependence of the lattice correlator 
G(t,T). Looking at Eq. (JXTJ) it is easy to see that the 
temperature dependence of G(r, T) comes from the tem- 
perature dependence of the spectral function o{uj,T) and 
the temperature dependence of the kernel K(t, ui, T). To 
separate out the trivial temperature dependence due to 
K(t,uj,T) following Ref. [17| we calculate the recon- 
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FIG. 8: The ratio G/Grecon for the pseudo-scalar channel on 
coarse £ = 2 lattices. 



FIG. 9: The ratio G/G reC on for the pseudo-scalar channel for 
the two finer £ = 4 lattices. 



structed correlator 

/•OO 

Grecon(T,T)= / dcja(u!, T = 0)K(t , cj, T). (35) 
Jo 

If the spectral function does not change with increas- 
ing temperature we expect G(r, 7 1 )/G recon (r, T) = 1. In 
this section we are going to study the temperature de- 
pendence of this ratio for different channels at different 
lattice spacings. To fix the temperature scale we use the 
(3 dependence of the r given by Eq. [3D] and the value 
roT c = 0.7498(50) [5(| for the transition temperature T c . 
The parameters of our finite temperature simulations are 
given in Tabic iLLTl 

A. The pseudo-scalar correlators 

First let us examine the temperature dependence of 
the pseudo-scalar correlators. In Fig. [8] we show our 
numerical results for G/G recon on coarse lattices with 
£ = 2. The figure shows very little temperature depen- 
dence of the correlators till temperatures 1.2T C . Calcu- 
lations at smaller lattice spacings enable us to consider 
higher temperatures. In Fig. [9] we show the tempera- 
ture dependence of G/G recon on our £ = 4 lattices. We 
see very little change in the pseudo-scalar correlator till 
temperatures as high as 1.5T C . Medium modifications of 
the correlator slowly turn on as we increase the tempera- 
ture above this value. From the figures it is clear that the 
temperature dependence of the correlators is not affected 
significantly by the finite lattice spacing. The very small 
temperature dependence of the pseudo-scalar correlator 
suggests that the corresponding ground state rj c (lS) sur- 
vives till temperatures as high a 1.5T C . The temperature 
dependence of the correlator found in this study is simi- 
lar to that of Ref. [l7| where isotropic lattices with very 
small lattice spacings, a -1 = 4.86, 9.72 GeV have been 
used. We find a somewhat stronger temperature depen- 
dence of G/G recon than in Ref. [l7j . In particular, at 



1.5T C we see small, but statistically significant deviations 
of G /G recon from unity. At higher temperatures the de- 
viation of G/ G recon from unity become slightly larger 
than those found in Ref. [13] ■ This is possibly due to the 
fact that cutoff effects are more important at higher tem- 
peratures. Thus despite similarities of the temperature 
dependence of the pseudo-scalar correlator to findings of 
Ref. we see quantitative differences. One should 

note, however, statistical errors and systematic uncer- 
tainties are larger in the analysis presented in Ref. [l7T | 
than in this calculation. The ratio G/G recon starts to 
depend more strongly on the temperature around 2T C . 
This may suggest some quantitative differences in the 
properties of the lowest state at this temperature. 



B. The P-wave correlators 

In this subsection we are going to discuss the temper- 
ature dependence of the scalar, axial-vector and tensor 
correlators corresponding to P-states. The numerical re- 
sults for the scalar correlator on our £ = 2 coarse lattices 
are shown in Fig. 1101 As one can see the correlator 
is temperature independent below T c and strongly en- 
hanced above T c . The magnitude of the enhancement 
is largest on the coarsest lattice and decreases with de- 
creasing lattice spacing. The numerical results on fine 
lattices are shown in Fig. QT] We see some differences 
in G /G recon calculated at (3 = 6.1 and (3 — 6.5. Thus 
the cutoff dependence of G/G recon is larger in the scalar 
channel than in the pseudo-scalar one. For (3 = 6.1 and 
£ = 4 we also did calculations on 24 3 x 24 lattice to 
check finite volume effects. The corresponding results 
are shown in Fig. [TT] indicating that the finite volume ef- 
fects are small. On the finest lattice the enhancement of 
the scalar correlator is very similar to that found in cal- 
culations done on isotropic lattices [l3], but small quan- 
titative differences can be identified. 

In Figs. [T2] and [T3] we show the temperature depen- 
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TABLE III: Simulation parameters for charmonium at finite temperature. We assumed roT c 
T c = 295(2) MeV. 
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FIG. 10: The ratio G/Grecon for the scalar channel on coarse 
£ = 2 lattices. 
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£ = 4 lattices. 
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FIG. 13: The ratio G/G re con for the tensor channel for £ 
lattices. 



dence of the axial-vector and tensor correlators respec- 
tively for £ = 4. Qualitatively their behavior is very 
similar to the scalar correlator but the enhancement over 
the zero temperature result is larger. The results for the 
axial-vector correlators again are very similar to those 
published in Ref. [l7|. The difference in G/G recon cal- 



culated at j3 — 6.1 and j3 = 6.5 are smaller than in the 
scalar channel. 

The large increase in the scalar, axial- vector and ten- 
sor correlators indicates strong modification of the corre- 
sponding spectral function and, possibly the dissolution 
of IP charmonia states. 
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C. The vector correlator 

The numerical results on the vector correlator are 
shown in Fig. [TH for £ = 4. As one can see from the 
figures the temperature dependence of G/G rec is differ- 
ent from the pseudo-scalar case and this ratio is larger 
than unity for all lattice spacings. Though this appears 
to be surprising at first glance it can be explained by 
the presence of the diffusion pole in the vector spectral 
function [ll|, [52J. Since the vector current is conserved 
there is a low energy contribution in the vector spectral 
function at finite temperature 



cr. 



(HI 



(36) 
(37) 



where (v 2 h ) is the average thermal velocity of the heavy 
quark. For free streaming heavy quarks at temperature 
T we have (v 2 h ) — T/M, with M being the heavy quark 
mass. Furthermore, Xs(T) is the quark number suscep- 
tibility corresponding to the heavy quark, which in the 
non-interacting case is given by: 



Xs(T) = 



T 



d 3 p 



1 



(2^) 3 exp((Vp 2 +Af 2 )/T) + 1 ' 



(38) 



The low energy parts of a^ and correspond to charge 
fluctuations and heavy quark diffusion, respectively. The 
smeared delta function, S v (uj), in Eq. (|36J) carries infor- 
mation about the heavy quark diffusion constant D [52| . 
Using effective Langevin theory [52j it can be shown that 



1 



V 



7T io 1 + r\ A 



V 



T 
AID' 



(39) 



Therefore a very precise calculation of the vector corre- 
lator at finite temperature can provide some information 
about heavy quark transport in Quark Gluon Plasma. 
This would require a method which is capable to iso- 
late reliably the low energy contribution given to the 
Euclidean correlator which is given by Eq. (|36|) . In 
the present paper we attempted to do so by subtracting 
from the finite temperature spectral functions the corre- 
sponding high energy part. We have assumed that this 
high energy part can be approximated reasonably well by 
the zero temperature spectral function. As we have seen 
in the previous subsection the pseudo-scalar correlators 
show little temperature dependence up to temperatures 
as high as 1.5T C . This suggests that to some extent the 
spectral function up to this temperature can be approx- 
imated by the zero temperature spectral function. Since 
the lowest lying states in the pseudo-scalar and vector 
channel are the same (up to the small hyperfinc splitting) 
we may expect that the same is true for the vector chan- 
nel, and that the only difference is the transport peak at 
u) ~ 0. Therefore we expect that G v — G^ econ should give 
an estimate for the low energy contribution coming from 
Eq. (|3"6"|) . It is easy to see that for the temporal compo- 
nent of the vector correlator we have Goo = —Txs{T) ( 
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FIG. 14: The ratio G/G r 
two finer £ = 4 lattices. 



for the vector channel for the 



c.f. Eq. (|37jl ). This is a consequence of charge conserva- 
tion. In Fig. [15] we show -(G v (t) - G^ econ (r))/G o as 
function of t. Note that this quantity does not depend 
on the renormalization constant. Since the width of the 
smeared delta function in Eq. (|36p is small the above 
quantity should be approximately independent of r and 
give the averaged thermal velocity (v 2 h ). The data in Fig. 
U5l do not agree completely with these expectations, devi- 
ations from this behavior are seen at short distances. The 
deviations are significant for 1.5T C . The problem is that 
with increasing temperature the high energy part of the 
spectral function starts to deviate from the zero temper- 
ature result, e.g. the properties of the lowest state could 
be slightly modified. The difference G y (r) - Gjf eco „(r) 
can be quite sensitive to these small changes, especially 
as the temperature increases. 

It has been noticed already in the previous subsec- 
tion that at 1.5T C we see small but significant tem- 
perature modifications of the ratio G/G recon for the 
pseudo-scalar channel. Therefore at this temperature 
we attempted to take into account possible modifica- 
tions of the high energy part of the spectral function. 
We do this by lowering the amplitude of the 1st peak 
by 7% in the zero temperature spectral function. In 
the pseudo-scalar channel this ensures that G/G recon 
stays around unity. Using this modified spectral func- 
tion in the vector channel we construct the quantity 
— (G v (t) — G y co „(T))/Goo again. This is also shown 
in Fig. [T3] as open symbols. We see that as a result of 
this correction ~(G v (t) — G,J4 co „(t))/Goo shows a nice 
plateau and also increases with the temperature as ex- 
pected. For the averaged thermal velocity squared we 
estimate (v 2 h ) ~ 0.13 for l.2T c and {v 2 th ) ~ 0.18 for 
1.5T C . A non-zero diffusion coefficient r\ would produce 
a small curvature in —(G v (t) — G^ co „(r))/Goo, clearly 
within present systematic uncertainties we cannot make 
any statement on the value of rj. 

The quantity (G — G recon ) / Goo gives an estimate of 
T/M and is shown in Fig. [T5J As one can see from the 
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FIG. 16: The dependence of the reconstructed pseudo-scalar 
spectral function on the maximal temporal extent for f3 = 6.5. 
In the analysis the default model m(u) = 1 has been used. 



figure it has the expected magnitude. 



VII. CHARMONIUM SPECTRAL FUNCTIONS 
AT FINITE TEMPERATURE 

The study of the temperature dependence of the char- 
monium correlators presented in the previous section pro- 
vides some useful insight into the properties of the spec- 
tral functions in the deconfined phase. In particular it 
suggests the survival of the IS state till temperatures as 
high as 1.5T C and significant modification or dissolution 
of the IP charmonium states. To get more detailed infor- 
mation about charmonium spectral functions we would 
like to use MEM. We attempted to reconstruct the char- 
monium spectral functions at finite temperature on our 
two finest lattices. Below we are going to present our 
results for charmonium spectral functions at finite tem- 
perature for different channels. 



A. Pseudo-scalar and vector spectral functions at 
finite temperature 

In section [V] we have seen that using MEM we can 
reconstruct well the main features of the spectral func- 
tion, in particular the ground state properties. At fi- 
nite temperature the situation becomes worse because 
the temporal extent is decreasing. The maximal time 
separation is r max = 1/(2T). As a consequence it is no 
longer possible to isolate the ground state well. Also the 
number of available data points becomes smaller. While 
the later limitation can be overcome by using finer and 
finer lattice spacing in time direction the former limita- 
tion is always present. Therefore we should investigate 
systematic effects due to limited extent of the temporal 



direction. It appears that the pseudo-scalar channel is 
the most suitable case for this investigation as at zero 
temperature it is well under control and there is no con- 
tribution from heavy quark transport. To estimate the 
effect of limited temporal extent we have calculated the 
spectral functions at zero temperature considering only 
the first Ndata time-slices in the analysis for (3 = 6.5, 
£ = 4. The result of this calculation is shown in Fig. [TBI 
where we considered Ndata — 80, 40, 20 and 16. The 
last two values correspond to our finite temperature lat- 
tices at this value of (3. In this case we see the first peak 
quite clearly when a — a max - For (3 = 6.1 as well as 
for (3 = 6.5 and T > 1.5T C the 1st peak is only clearly 
visible when we consider values for a which are smaller 
than a max - As one can see from the figure already for 
Ndata — 40 and T max — 0.56 fm the second peak corre- 
sponding to radial excitation is no longer visible and the 
first peak becomes significantly broader. The position of 
the first peak, however, is unchanged. As the number of 
data points is further decreased (Ndata — 20, 16) we see 
further broadening of the first peak and a small shift of 
the peak position to higher energies. These systematic 
effects should be taken into account when analyzing the 
spectral functions at finite temperature. Therefore when 
studying the spectral functions at finite temperature we 
always compare with the zero temperature spectral func- 
tions reconstructed with the same number of data points 
and T max as available at that temperature. 

In Fig. [T7] the spectral functions at different tempera- 
tures are shown together with the zero temperature spec- 
tral functions. We used m(u>) — 0.01 as a default model 
for T = 1.2T C and T = 1.5T C . For T = 2.0T C we used 
m(u) = 1 since the use of m(uo) = 0.01 resulted in mul- 
tiple maxima for P[a\DH]. The figure shows that the 
pseudo-scalar spectral function is not modified till 1.5T C 
within the errors of our calculations. This is consistent 
with the conclusions of Ref. [ljl [13] ■ One should note, 
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FIG. 17: The pseudo-scalar spectral function for /3 — 6.5 and 
Nt — 40, 32,20 corresponding to temperatures 1.2T C , 1.5T C 
and 2.4T C . In the analysis the default model m(ui) — 0.01 
has been used for 1.2T C and 1.5T C , while for 2.0Tc we used 
m(w) = 1. 



however, that it is difficult to make any conclusive state- 
ment based on the shape of the spectral functions as 
this was done in the above mentioned works. The de- 
pendence of the reconstructed spectral functions on the 
default model m(o>) is much stronger at finite tempera- 
ture. We have reconstructed the spectral functions using 



different type of default models. For all temperatures 
T < 1.5T C the difference between the finite temperature 
spectral function and the zero temperature one is very 
small compared to the statistical errors for all default 
model considered here. We discuss this further in the 
Appendix. In particular we used the default model con- 
structed from the high energy part of the lattice spectral 
functions calculated at zero temperature as this was done 
in Ref. [l7|. The idea is that at sufficiently high energy 
the spectral function is dominated by the continuum and 
is temperature independent and is suitable to provide the 
prior knowledge, i.e. the default model. This behavior 
of m(uj) is then matched to a simple uj 2 dependence [l7l ]. 
With this default model we have calculated the spectral 
function also at /3 = 6.1 on 16 3 x 26 and 16 3 x 24 lattices 
corresponding to 1.07T C and 1.16T C . Also here we have 
found very little temperature dependence of the spectral 
functions. The survival of the IS charmonium state was 
also confirmed by the analysis of the correlation functions 
with different spatial boundary condition (53j. 

At sufficiently high temperatures IS" charmonium 
states should dissolve and thus we expect significant 
changes in the spectral functions. The analysis of the 
pseudo-scalar correlators suggests that this may happen 
around 2T C which was also the conclusion of Ref. jig . 
The spectral function at T = 2.0T C shown in Fig. [T7] 
shows significant deviation from the zero temperature 
spectral function. However, this conclusions could be 
premature, since for a different default model, namely 
the free lattice spectral function, we see almost no tem- 
perature dependence of the spectral function also at 2T C . 
Using the high energy part of the zero temperature spec- 
tral function as the default model as discussed above we 
arrive at a similar conclusion. Thus further studies are 
needed to establish at which temperature the 15 char- 
monia states dissolve. 

We also calculated the spectral function in the vector 
channel. The results are shown in Fig. [18] for the de- 
fault model m(uj) = 0.01 . As this was already discussed 
in the previous section the basic difference between the 
pseudo-scalar and vector spectral functions at finite tem- 
perature is the presence of the transport peak at uj ~ 0. 
The difference of the temperature dependence of the vec- 
tor and pseudo-scalar correlators is consistent with this 
assumption. The vector spectral function reconstructed 
with MEM shows no evidence of the transport peak at 
uj ~ 0. On the other hand the spectral function at 1.2T C 
differs from the zero temperature spectral function, in 
particular the first peak is shifted to smaller uj values. 
We believe this is a problem of the MEM analysis which 
cannot resolve the peak at w ~ but instead mimics its 
effect by shifting the J ftp peak to smaller uj. Also at 2.4T C 
the spectral function extends to smaller uj values than in 
the pseudo-scalar correlator which again indicates some 
structure at a; ~ 0. We analyzed the vector spectral 
functions using other choices for the default model and 
always found that the spectral functions at finite temper- 
ature differs from the zero temperature spectral functions 
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FIG. 18: The vector spectral function for (3 — 6.5 and Nt = 
40, 20 corresponding to temperatures 1.2T C and 2.4T C . In the 
analysis the default model m(uj) — 0.01 has been used. 



and extend to significantly smaller oj values. 
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FIG. 19: The scalar spectral function for f) = 6.1 at T = 
1.16T C and at zero temperature reconstructed using Ndata = 
12. At finite temperature two default models m(ui) — 0.01 
and m(uj) = 0.038lj 2 have been used. 



fled compared to the zero temperature spectral function; 
its shape at small u>, however, strongly depends on the 
default model. We also analyzed the spectral function for 
these data with the default model constructed from the 
high energy part of the zero temperature lattice spectral 
functions as described in the previous subsection. The 
spectral functions obtained from this analysis are similar 
to those shown in Fig. [19] (see Appendix). 

From the analysis of the spectral functions we conclude 
that the IP states are dissociated in the deconfined phase 
at temperatures T ~ 1.1 — 1.2T C . This is consistent with 
the analysis of the correlators presented in the previous 
section as well as with the conclusions of Ref . [17| ■ On 
the other hand with present level of statistical accuracy 
we cannot make precise statements about the form of the 
P-wave spectral function at high temperatures. 



B. The P-wave spectral functions 

We have studied the spectral functions in the scalar, 
axial- vector and tensor channels. Already at zero tem- 
perature the calculations of the P-wave spectral functions 
turned out to be much more difficult than for the S-wave 
spectral functions. This is even more the case at finite 
temperature where the number of data points is limited. 
For most of the choices of the default model we have 
found that P[a\DH] has multiple maxima for all of our 
data sets. This is a sign of insufficient statistics. The only 
exception is the data set corresponding to 16 3 x 24 lattice 
with f3 = 6.1, i.e. T = 1.16T C . The spectral function in 
the scalar channel for j3 = 6.1 at T = 1.1 QT C is shown 
in Fig. [HI] and compared with spectral functions at zero 
temperature reconstructed using Ndata — 12 data points. 
At finite temperature two default models, m(ui) — 0.038 
and m(uj) — 0.01 have been used. We see that the spec- 
tral function at finite temperature is significantly modi- 



VIII. CHARMONIUM CORRELATORS AND 
SPECTRAL FUNCTIONS AT FINITE MOMENTA 

So far we considered charmonia at zero spatial momen- 
tum, i.e. charmonia at rest in the heatbath's rest frame. 
It is certainly of interest to study the temperature de- 
pendence of correlators and spectral functions at non- 
zero spatial momentum. Such a study has been done us- 
ing isotropic lattices with lattice spacing a" 1 = 4.86GeV 
and 9.72 GeV [13, Hi]. It has been found that the pseudo- 
scalar correlators are enhanced compared to the zero tem- 
perature correlators for non-vanishing spatial momenta. 
We have calculated the finite momentum pseudo-scalar 
correlators at j3 = 6.1, £ = 4 for different temperatures. 
The results are shown in Fig. [501 We see that the temper- 
ature dependence of G/G recon at zero and finite momen- 
tum are different. At zero momentum this ratio decreases 
monotonically with increasing r and increasing temper- 
ature. At finite momentum this is not the case, we see 
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that the decrease in G/G recon is at least weaker, and at 
sufficiently large momenta we see even that this ratio in- 
creases with r. The differences in G/G re con calculated in 
this work and in Refs. [17l.l54| are present already at zero 
momentum and are presumably due to finite lattice spac- 
ing errors. Apart from this the momentum dependence 
of the pseudo-scalar correlators is similar to the findings 
of Refs. [13, HH ■ We also calculated the spectral function 
at finite temperature and found that the temperature de- 



pendence of the pseudo-scalar spectral functions for non- 
zero momenta is significant. This again is in qualitative 
agreement with the results based on isotropic lattice cal- 
culations. It would be interesting to see if the difference 
in the temperature dependence of the correlators at zero 
and finite spatial momenta is due to a contribution to 
the spectral functions below the light cone at finite tem- 
perature unm- 



ix. BOTTOMONIUM SPECTRAL FUNCTIONS 
AT ZERO TEMPERATURE 

The use of Fermilab formulation described in the pre- 
vious sections allows us to study also bottomonium for 
the same range of lattice spacings. Usually bottomonium 
is studied using lattice NRQCD (see e.g Ref. |55j). So 
far the only study of bottomonium within the relativistic 
framework was presented in Ref. [5l[ . The simulation pa- 
rameters for bottomonium are summarized in Table [TVl 
As before the lattice spacing is fixed by the Sommer scale 
ro, assuming vq = 0.5fm. The values of ro for f3 = 5.9 
and 6.1 are taken from Ref. [46j], while for (3 — 6.3 we 
use the extrapolation described in section HVl The bare 
quark mass, the bare velocity of light v t as well as the 
clover coefficients were taken from Ref. 5l|. Note that 



the lattice spacings determined from ro are about 20% 
smaller than in Ref. [5l| where it was determined from 
bottomonium X P\ — 15 mass splitting. As the result our 
estimates of the T mass are smaller than those in Ref. 
[5l( | . However, if we use the values of the lattice spacing 
quoted in Ref. [51] to calculate the physical masses we 
find good agreement. 

In our study we will stick to the value of the lat- 
tice spacing determined from ro as it agrees reasonably 
well with the lattice spacings extracted from the charmo- 
nium 1 P 1 — 15 splitting. This splitting is apparently less 
sensitive to the effect of quenching. The bottomonium 
1 P\ — 15 splitting, on the other hand, is more sensitive 
to quenching errors as well as more difficult to estimate 
precisely on the lattice. In Table Uvl we also show the T 
masses. 

Using MEM we analyzed the spectral functions in dif- 
ferent channels for all three lattice spacings. In Fig. [21] 
we show the spectral functions in the pseudo-scalar chan- 
nel. Since the physical quark mass is different at differ- 
ent lattice spacings we shifted the horizontal scale by 
the difference of the calculated T-mass and the corre- 
sponding experimental value. We can see that the first 
peak in the spectral function corresponds to the rjb{XS) 
state and its position is independent of the lattice spac- 
ing. The remaining details of the spectral functions are 
cut-off dependent and we cannot distinguish the excited 
states from the continuum. The position and the ampli- 
tude of the first peak in the spectral functions is in good 
agreement with the results of simple exponential fit. As 
in the charmonium case the maximal energy Lo max for 
which the spectral function is non-zero scales approxi- 
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TABLE IV: Simulation parameters for bottomonium at zero 
temperature 
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FIG. 21: The pseudo-scalar bottomonium spectral function 
at zero temperature for different lattice spacings. 



mately as a s 1 . Similar results have been obtained in the 
vector channel. 

The spectral function in the scalar channel is shown 
in Fig. [22] As it was the case for charmonium the 
correlators in this channel are more noisy than in the 
pseudo-scalar channel and as a result it is more difficult 
to reconstruct the spectral function. Nevertheless we are 
able to reconstruct the Xbo state which is the first peak 
in the spectral function. The peak position and the am- 
plitude are in reasonable agreement with the result of 
simple exponential fit. 
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FIG. 22: The scalar bottomonium spectral function at zero 
temperature for different lattice spacings. 



X. BOTTOMONIA AT FINITE TEMPERATURE 

Having calculated the bottomonium spectral function 
at zero temperature we can study the temperature depen- 
dence of bottomonium correlators to see medium modifi- 
cation of bottomonia properties. The parameters of the 
finite temperature simulations are summarized in Table 

m 

In Fig. [53] we show G/G recon for vector and pseudo- 
scalar channel at different lattice spacings. This ratio 
appears to be temperature independent and very close 
to unity up to quite high temperatures. This is con- 
sistent with the expectation that IS bottomonia are 
smaller than IS charmonia and thus less effected by the 
medium. They could survive till significantly higher tem- 
peratures. Compared to charmonium case the differ- 
ence between the pseudo-scalar and vector channels is 
smaller. This is also expected as the transport contribu- 
tion which is responsible for this difference is proportional 
to Xs(T) ~ exp(— m c ,b/T), and thus is much smaller for 
bottom quarks. 

The temperature dependence of the scalar correlator 
is shown in Fig. [2U Contrary to the pseudo-scalar 
and vector correlators it shows strong temperature de- 
pendence and G/G rec0 n is significantly larger than unity 
already at 1.1T C . Since the size and the binding energy 
of IP bottomonium states are comparable to that of the 
IS charmonium states one would expect a much smaller 
temperature dependence. The data on the other hand 
shows temperature dependence which is even larger than 
for IP charmonium states. This may suggest that the 
IP bottomonium states are strongly modified or maybe 
even dissolved at T ~ 1.1 — 1.2T C . Further investigation 
are needed to clarify what happens to IP bottomonium 
states above the deconfinement transition. Model studies 
of the P-wave bottomonium correlators suggest that its 
large change may be due to the modification of the con- 
tinuum part of the spectral functions and not to the dis- 
solution or strong modification of the ground state [Til ] . 
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TABLE V: Simulation parameters for bottomonium at finite temperature. We assumed roT c = 0.7498(50) corresponding to 
T c = 295(2) MeV. 



CD 

<3 



CD 

<3 



CD 

<3 



CD 

<3 



CD 

<3 



CD 

<3 




; ■ 

I t * » i i i ; ; * 



0.1 



0.2 
x[fm] 



0.3 



FIG. 23: The ratio G/G rec in the pseudo-scalar (top) and 
vector channels at different lattice spacings. 



We have tried to reconstruct the spectral functions at 
finite temperature for (3 = 6.3. For the pseudo-scalar 
channel the results are shown in Fig. [25] As in the char- 
monium case we compare the finite temperature spec- 
tral function with the zero temperature spectral func- 
tion obtained with the same number of data points and 
time interval. As expected the spectral function shows 
no temperature dependence within errors. On the other 
hand we could not reliably reconstruct the scalar spectral 
function at finite temperature, the probability P[cv|I?to] 
showed no maximum as the function of a. Presumably 
much more statistics is needed to get some information 
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FIG. 24: The ratio G/G rec in the scalar channel for different 
lattice spacings. 



about the scalar spectral function. 



XI. CONCLUSIONS 

In this paper we investigated quarkonium correlators 
at zero and finite temperature in quenched QCD using 
anisotropic lattices and the Fermilab formalism for the 
heavy quarks. In our investigations we used a variety of 
different lattice spacings to control cut-off effects in the 
quarkonium correlators and spectral functions. Using the 
Maximum Entropy Method we have calculated the spec- 
tral functions both at zero and finite temperatures. This 
was done using a newly developed algorithm for the MEM 
analysis which does not rely on singular value decompo- 
sition. It appeared to be more stable than the Bryan al- 
gorithm used so far in the analysis of the meson spectral 
functions. Therefore we could for the first time reliably 
calculate the charmonium spectral functions at zero tem- 
perature identifying the main structures : ground state, 
excited states and the continuum. Using the zero tem- 
perature results as the base line we investigated the tem- 
perature dependence of charmonium correlators at finite 
temperature at different lattice spacings and showed that 
it is qualitatively the same for all lattice spacings. We 
studied the charmonium spectral functions at finite tem- 
perature paying attention to systematic effects arising 
from the fact that the time interval at finite tempera- 
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FIG. 25: The pseudo-scalar bottomonium spectral function 
at finite temperature. 



ture as well as the number of data points decreases with 
increasing temperature. We found that the spectral func- 
tions in the pseudo-scalar channel do not change up to 
temperatures 1.5T C within systematic and statistical er- 
rors of our calculations. We have also learned that it is 
difficult to make statements about the survival of char- 
monium states based solely on the shape of the spectral 
function as it is distorted by the limited time extent and 
strongly depends on the default model. On the other 
hand the small temperature dependence of the pseudo- 
scalar spectral function up to 1.5T C holds for all default 
models. The correlation functions in all channels corre- 
sponding to IP states show very different behavior. To- 
gether with the results for the scalar spectral functions 
this suggests the melting of the IP charmonium states 
at temperature T = 1.1 — l.2T c . In general, most of our 
findings in charmonium sector is in reasonable agreement 
with results known from studies on fine isotropic lattices 
[l7| . We also identify the transport contribution in the 
vector correlators for the first time. 

Our investigations in the bottomonium sector are less 
detailed. We were able to reconstruct the bottomonium 
spectral functions at zero temperature and show that the 
first peak is not affected by lattice artifacts. Nonetheless 
the spectral functions are less accurate than for charmo- 
nium. We were not able to reconstruct reliably the bot- 
tomonium spectral functions at finite temperature. The 
reason is that because of the large quark mass the corre- 
lators decay fast and the signal is poor. The temperature 
dependence of the correlators suggests that the IS bot- 
tomonium state can survive till much higher temperature 
than the charmonium states. On the other hand we see, 
unexpectedly, drastic changes in the IP bottomonium 
correlators right after deconfinement transition. To clar- 
ify the fate of the IP bottomonium states further stud- 
ies are clearly needed. They presumably should rely on 
NRQCD where the large bottom quark mass is integrated 
out. 

The use of anisotropic lattices and Fermilab approach 



reduces discretization effects and allows to study the 
problem of quarkonium dissolution on coarser lattices. 
In fact the analysis of the temperature dependence of 
the correlation functions shows that they can be reliably 
studied already for lattice spacings a^ 1 = 2.0GcV. Al- 
though the temperature dependence of the correlators 
turn out to be the same for all lattice spacings studied in 
this paper, we clearly see quantitative differences. This is 
important when we want to extend the studies of quarko- 
nium spectral functions to full QCD, where the lat- 
tice spacing is main limiting factor. Preliminary results 
from calculations of charmonium spectral functions in 
full QCD have been reported in Ref. 56]. These results 
are in qualitative agreement with findings in quenched 
QCD. However, no quantitative analysis has been per- 
formed so far to estimate the effect of see quarks. To get 
sufficiently close to the continuum limit in the analysis of 
the correlation function the spatial lattice spacing should 
presumably be smaller than (4GeV) _1 , i.e. it should be 
in the range used in studies with isotropic lattices (l7j . 
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Appendix 

In this appendix we present some technical details on 
the calculations of the quarkonium spectral functions. To 
study the reliability of MEM at zero temperature com- 
pare the masses and amplitudes of different states ob- 
tained from MEM and simple 2-exponential fit. This 
comparison is shown in Table IVII We also tried to esti- 
mate the systematic errors in the fit by allowing higher 
excited states (3 exponential fit). We show this system- 
atic uncertainty in the square brackets. One can see that 
we have a good agreement for the ground state in pseudo- 
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FIG. 26: The pseudo-scalar spectral function for different val- 
ues of a. The spectral functions were calculated on 16 3 x 96 
lattice at = 6.1, £ = 4. 
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FIG. 27: The conditional probability P[a|Dm] for different 
algorithms and different t max corresponding to the pseudo- 
scalar spectral function calculated on 16 3 x 96 lattice at = 
6.1, £ = 4. 



scalar channel (r] c ), the small disagreement in the values 
of the meson masses (typically appearing in the 4th digit) 
is the consequence of the discretization in lo. We have 
checked that for sufficiently fine binning of the spectral 
function this remaining error could be removed. For the 
scalar channel we find agreement between the results of 
the fit and MEM within estimated errors, with the only 
exception being = 6.1 where we find small deviation 
for the amplitudes which are outside the errors. In the 
pseudo-scalar channel we find also an agreement for the 
masses and amplitudes of the excited states within the 
estimated errors. 

To compare the spectral functions at different lattice 
spacings one needs the renormalization constants of dif- 
ferent quark bilinears. For the lattice fermions used in the 
present paper these have not been calculated. Therefore 
we introduced ad-hoc renormalization constants which 
are shown in Table IVII1 

The real parameter a enters into MEM through Eq. 
(jT2J . Therefore we have studied the dependence of the 
spectral functions on a. For the pseudo-scalar channel 
at zero temperature this is shown in Fig. [26j To get 
rid of the alpha dependence of the spectral function one 
calculates the conditional probability P[a\DH]. For suf- 
ficiently good quality Monte-Carlo data this function has 
a unique maximum at some alpha. This is shown in Fig. 
1271 for the new algorithm used in this paper as well as 
for the Bryan algorithm. The spectral functions calcu- 
lated with the new algorithm and the Bryan algorithm 
are shown in Fig. [28] 

At finite temperature the dependence of the spec- 
tral functions on the default model becomes important. 
Therefore we investigated the default model dependence 
of the spectral function. First we used different func- 
tional forms for the default model, m(uj) ~ uj 2 motivated 
by the free spectral functions of the continuum QCD and 
m (w) = const. We also considered the free lattice spec- 
tral function calculated in Ref. |24| as a default model. 
In Fig. [25] we show the pseudo-scalar spectral function 
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FIG. 28: The pseudo-scalar spectral function calculated on 
16 3 x 96 lattice at = 6.1, £ = 4 using the algorithm de- 
scribed in this paper and the Bryan algorithm for a = a max . 
Also shown in the figure is the spectral function obtained by 
integrating over a with the weight P[a|Dm]. See the main 
text for further details. 



at 1.2T C calculated with different default models. The 
figures shows that the default model dependence of the 
spectral function is significant, in particular the first peak 
is not present for all default models. Following Ref. [l7| 
we have used the high energy part of the zero tempera- 
ture lattice spectral function and matched it to the u> 2 be- 
havior at some energy uj t h above the energy region where 
no individual resonances can be seen. In the pseudo- 
scalar channel at = 6.5 we have chosen cu t h — 7GeV, for 
— 6.1 we have chosen ui t h — 5GeV in the pseudo-scalar 
channel and u t h — 5.5 GeV in the scalar channel. In Fig. 
[30] we show the scalar spectral function at zero temper- 
ature with Ndata — 12 and at T = 1.16T C for this type 
of the default model. The scalar spectral function shows 
significant change with the temperature in this case too. 
We analyzed the pseudo-scalar spectral function with this 
type of the default model and the result of this analysis 
is shown in Fig. [3U For this choice of the default model 
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0.281(8) 


0.285(2) 



TABLE VI: The zero temperature charmonium masses and amplitudes obtained from MEM and compared to 2-exponential fit 
with jackknife errors. We show the systematic errors of the fit procedure in the square brackets. 



0, c 


6.1, 2 


6.1, 4 


6.5, 4 


Zps 


0.374 


0.447 


0.447 


Zvc 


0.343 


0.447 


0.426 


Zsc 


0.534 


0.791 


0.620 


Zax 


0.592 


0.949 


0.707 



TABLE VII: The ad-hoc renormalization constants for char- 
monium used in Figs. 2, 3, 5-7, 15-19 in different channels. 
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free lattice 
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"a 
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FIG. 29: The default model dependence of the pseudo-scalar 
spectral function at 1.2T C calculated on 24 3 x 40 lattice with 
13 = 6.5. 




significant temperature dependence can be seen only for 
T = 2AT C . In the figure we also show the errors where 
appropriate. We find that the statistical errors are much 
smaller for this default model. 




10 12 14 
<o[GeV] 

FIG. 30: The scalar spectral function at 1.16T C and at zero 
temperature reconstructed using Tmax = 12 calculated with 
the default model coming from the high energy part of the 
zero temperature spectral function (see text). Also shown 
the zero temperature spectral function calculated using all 
data points and m(w) = 1. 
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